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EFFICIENT  CALCULATION  OF  DIRECTIVITY  INDICES 
FOR  CERTAIN  THREE-DIMENSIONAL  ARRAYS 

INTRODUCTION 

Directivity  Index  (DI)  may  be  defined  as  a  measure  of  the  improvement  in  the  signal-to- 
noise  power  ratio  (S/N)  that  a  beamformed  array  provides  in  an  ideal  isotropic  noise  field  with  a 
perfectly  correlated  signal,  relative  to  an  omni-directional  array  element  in  the  free-field,  i.e., 

0/  ^  10  log  OF  =  10  Lost  - ■  (« 

^  '  'free- field  omni-element 

For  an  array  with  arbitrary  geometry  and  element  locations,  the  beamformed  array  single- 
frequency  amplitude  response^  may  be  written  as 

A(0,(f))  =  Xwj  gj(0,(l))  exp(i(k-?g)  •rj) ,  (2) 

t=l 

where  A(0,(|))  is  the  array's  angular  response  to  an  arrival  firom  azimuthal  angle  0  and  polar  angle 
(j),  T  is  the  total  number  of  array  elements,  w^  and  g^  are,  respectively,  the  amplitude  shading 
and  element  angular  sensitivity  of  the  t-th  array  element,  k  is  the  wavevector  corresponding  to 
an  arriving  acoustic  planewave,  and  Ics  is  the  wavevector  to  which  the  array  is  steered.  The 
corresponding  steering  angles  are  0^  ,<j)g .  The  t-th  array  element  position  is  defined  by  the 
coordinate  vector 

A  spatial  coordinate  system  (with  corresponding  wavevectors)  for  an  array  oriented  in  the 
xz-plane  is  illustrated  in  figure  1,  along  with  the  definition  of  spherical  angles  0,(1).  The 
programs  supplied  in  the  appendices  of  this  report,  to  calculate  directivity,  use  the  geometry  of 
figure  1  to  define  input  parameters.  For  the  programs  planar-xz-equal.f  and  planar-xz-grid.f,  a 
rectangular  array  in  the  xz-plane  is  assumed.  Program  planar-xz-grid.f  allows  for  arbitrary 
element  grid  spacing,  that  is,  dx  or  dz  may  vary  from  element-to-element  The  program 
voIumetric-grid.f  assumes  an  array  configuration  which  allows  for  cylindrical  curvature  about 
the  x-axis.  At  each  fixed  array  row  x^ ,  elements  are  located  at  the  same  positions  {yn,  Zn);  if  all 
yn  are  equal,  the  array  becomes  a  plane  parallel  to  the  xz-plane. 
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Figure  1.  Array  Coordinate  System  and  Spherical  Angles  Denoting 
Incident  Planewave  Arrival  Direction  6,<|> 

With  the  aid  of  equations  (1)  and  (2),  the  components  which  comprise  the  calculation  of 
the  directivity  index  can  be  written  as 


5^  .  =S, 

ff—omiu 

where  S  is  defined  as  the  given  signal  power  level,  while 

2x  K 

N,^  =  \\Nsm(4,)d<pde, 

0  0 

where  N  is  defined  as  the  isotropic  noise  level,  while 


(3a) 


(3b) 

(3c) 
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lizit 

Narray  =  j  j  N\A(e,<l>f  Sin(<{>)d(l)dd  .  (3d) 

0  0 

Noting  that  equation  (3c),  is  simply  47cN,  and  canceling  the  signal  and  isotropic  noise 
levels,  equation  (1)  reduces  to 


DI  =  10  Log 


2%  It 

J  J|A(e,(|))psin((l))d(|)de 

.  0  0 


2 

where  |A(0g,0g)  is  the  array  power  response  at  steering  angles 


(4) 


Equation  (3d)  indicates  an  integral  with  limits  that  account  for  all  acoustic  arrival  angles. 
Hence,  the  directivity  expression  above  assumes  an  unbaffled,  or  fiee-field,  array  configuration. 
In  a  baffle  configuration,  would  be  reduced,  resulting  in  an  increased  directivity. 

However,  the  programs  supplied  here  include  directional  (such  as  cosine  to  a  power)  element 
sensitivities,  gj(0,<l))  in  (2),  which  inherently  provide  the  gains.realized  from  a  baffle.  There  is 
no  need  to  correct  the  directivity  values  yielded  by  the  programs  to  account  for  baffling  of 
isotropic  free-field  noise. 

For  a  line  array  with  constant  inter-element  spacing,  the  double  integral  in  equation  (4) 
may  be  evaluated  exactly  at  the  frequency  for  which  spacing  dx  =  X/2,  regardless  of  the  steering 
direction.  This  leads  to 


DI  =  \0\o%\ 


f  T 

Vt=l  J 


t=l 


(5) 


Accurate  directivity  predictions  for  arbitrary  array  configurations,  at  general  frequencies 
of  interest,  require  numerically  integrating  the  double  integral  given  in  equation  (4)  and 
calculating  the  maximum  array  power  response,  ‘  should  be  noted  that  the 

maximum  array  response  may  not  occur  precisely  at  the  steering  angles  (^g.^^)-  Array 
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curvature,  unequal  array  element  spacing,  element  sensitivity,  and  shading  weights  act  to  shift 
the  maximum  response  axis  from  to  The  programs  supplied  here  have  global 

and  local  search  procedures  which  obtain  the  peak  power  response  and  the  corresponding  angles 
(d  ,0  ). 

The  remainder  of  this  report  wiU  discuss  accurate  and  efficient  numerical  integration 
methods  for  evaluating  the  double  integral  of  equation  (4).  The  focus  is  to  reduce  the  number  of 
evaluations  required  of  the  integrand,  in  order  to  achieve  rapid  convergence  and  accuracy  with  a 
minimum  of  computational  effort. 

For  very  small  partitions,  the  Riemann  sum  will  provide  a  reasonable  approximation  of 
the  definite  integrals  in  equation  (4).  Indeed,  if  an  exceedingly  large  number  of  narrow 
rectangles  (or  augmented  partitions)  are  used  in  the  Riemann  sum,  the  approximation  would  be 
quite  accurate.  This  implies,  however,  an  enormous  amount  of  function  evaluations,  each  of 
which  is  computationally  time-consuming.  For  a  large  number  of  array  elements,  on  the  order  of 
thousands,  such  a  simplistic  approach  would  be  unfeasible.  Vector  and  matrix  operations  used  in 
high-performance  computation  software,  such  as  MATLAB®,  do  provide  quick  function 
evaluations  and  have  built-in  numerical  integration  routines.  However,  these  computational 
benefits  can  be  reduced  by  the  amount  of  RAM  available  on  a  given  computer  system. 

Swapping  memory  space  between  a  system’s  hard  drive  and  RAM  can  be  inefficient;  resulting  in 
lengthy  integration  times.  The  goal  then  is  to  use  an  integration  method  which  yields  a  highly 
accurate  integral  approximation  with  a  minimum  number  of  function  evaluations. 

rx 

Due  to  the  periodicity  and  analyticity  of  the  integrand  |A(0,<t))|  in  azimuthal  angle  0,  an 
elementary  closed-type  Trapezoidal  formula  was  used  for  evaluating  the  outer  0  integral  in 
equation  (4)  over  a  full  period  (0,2n).  The  timer  integral,  over  polar  angle  <j),  uses  a  Gauss- 
Legendre  quadrature  formula  of  order  m,  where  m  can  vary  over  the  twelve  values  (16, 24,  32, 
48,  64, 96, 128, 192, 256,  384, 512, 768).  The  Gauss-Legendre  method  is  generally  significantly 
more  accurate  than  an  equal-interval  formula,  such  as  the  Trapezoidal  rule,  for  a  given  number  of 
function  evaluations.  However,  Gauss  quadrature  was  deliberately  not  chosen  for  the  outer  0- 
integral  due  to  the  particular  integration  limits  and  the  integrand  periodicity  and  analyticity  in  0. 
The  Trapezoidal  rule  is  extremely  efficient  for  integration  of  an  analytic  periodic  integrand,  when 
conducted  over  a  full  period.^  Generally,  efficient  numerical  evaluation  of  multi-dimensional 
integrals  requires  examination  of  the  nature  of  the  integrand  and  the  integration  limits;  rarely  is 
one  procedure  optimum  for  all  integrals. 
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FACTORIZATION  OF  ARRAY  RESPONSE 


The  array  amplitude  response  A(9,<1))  was  given  in  equation  (2)  for  the  general  three- 
dimensional  array  with  arbitrary  element  locations,  weights,  and  element  responses,  where  T  is 
the  total  number  of  elements  in  the  array.  In  the  general  case,  the  sum  of  T  terms  can  be 
evaluated  only  by  calculating  every  individual  complex  term  and  summing  up  those  T  quantities. 
This  is  a  particularly  time-consuming  task,  especially  when  T  is  of  the  order  of  10,000  or  more. 
When  coupled  with  the  fact  that  the  directivity  index  requires  computation  of  a  double  integral, 
namely  equation  (4),  which  must  be  repeated  for  each  different  frequency  and/or  steering  angle, 
the  computational  burden  becomes  excessive.  Furthermore,  as  the  number  T  of  elements 
increases,  the  array  power  response  |A(0,<l))j  becomes  even  sharper  in  angles  thereby 
requiring  still  finer  evaluation  of  the  integrand  in  equation  (4),  in  order  to  retain  accuracy  in  the 
final  DI  calculation. 

In  this  section,  we  will  significantly  reduce  this  computational  burden  for  a  particular 
class  of  arrays  and  weights,  by  factoring  the  amplitude  response  A(0,<|>)  into  a  product  of  two 
sums  which  have  far  fewer  total  terms  than  T.  For  example,  an  array  of  50  by  200  elements  will 
require  evaluation  of  50  +  200  =  250  terms  instead  of  50  *  200  =  10,000  terms,  a  savings  in 
execution  time  by  a  factor  of  40.  Additional  savings  in  time  are  achieved  by  precomputing  and 
storing  quantities  that  are  used  repeatedly  in  the  program  for  calculation  of  the  DI.  This  storage 
requirement  is  minimal,  even  for  the  very  large  arrays  of  interest  here. 


ELEMENT  LOCATIONS 

We  use  the  coordinate  system  illustrated  in  figure  1.  Receiving  elements  will  be  located 
only  on  L  planes  of  fixed  x,  namely,  at  {x^ }  for  1  ^  ^  <  L,  where  each  x^  is  arbitrary.  Thus, 
these  planes  need  not  be  equally  spaced  in  x. 

At  each  x-coordinate  x^ ,  N  elements  are  then  located  at  the  points  yn,Zn  for  1  <  n  <  N; 
these  pairs  {yn,Zn}  are  identical  for  all  x^ .  However,  each  location  yn.zn  is  arbitrary.  This 
generality  allows  for  the  array  to  take  on  a  cylindrical  shape,  for  example;  other  more  general 
three-dimensional  arrays  are  also  allowed.  In  the  special  case  where  the  {yn}  are  aU  equal,  this 
becomes  a  planar  array,  parallel  to  the  xz-plane;  however,  the  L  locations  {x^}  and  the  N 
locations  {zn}  are  still  arbitrary. 


5 


The  array  is  seen  to  contain  a  total  of  T  =  L  N  elements,  in  L  parallel  planes  of  N 
elements  each,  with  an  identical  repeated  (arbitrary)  configuration  in  each  x-plane.  This  is  the 
general  configuration  considered  here.  Any  three-dimensional  array  which  cannot  be  put  into 
this  form  will  not  be  covered  by  the  following  analysis  and  simplifications. 

Each  element  at  has  its  normal  in  the  yz-plane,  at  an  angle  \|/n  relative  to  the 

z-axis;  see  figure  2.  That  is,  the  element  lies  in  the  yz-plane  for  every  element  in  the  array; 
however,  the  element  orientation  angle  \|/n  can  vary  with  location  yn,Zn,  but  not  with  location  x^ . 
Thus,  the  direction  vector  of  the  n-th  element  boresight  is  0  i  H-  sin(\}rj^)  j  +  cos(\i/jj)  k  for  the 
element  at  x^,yj^,z^. 


Figure  2.  Element  Locations  at  a  Fixed  x^ 
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ELEMENT  RESPONSES 


The  cosine  of  the  angle  Bn(0,(}))  between  the  n-th  element  normal  and  a  general  arrival 
direction  0,({)  is  given  by  the  dot  product 

cos  (0,  (j))  =  [0  i  +  sin(\|/jj )  j  +  cos(\i/„ )  k]  • 

•  [cos(0)  sin((())  i  +  sin(0)  sin((j))  j  +  cos(<{))  k]  =  (6) 

=  sinCxjTjj)  sin(0)  sin(<t))  +  cosCxj/j^)  cos(<t)) . 

Notice  that  this  quantity  is  independent  of  i  and  ,  due  to  the  element  locations  and  boresight 
orientations. 

We  are  interested  here  in  relative  element  responses  of  the  form 

rn(0,<t))  =  [cosB^(0,(t))f  for|B^(0,(l))|^7c/2,  (7) 


where  d  is  a  power  under  our  control.  By  means  of  the  procedure  presented  in  Appendix  A,  the 
element  responses  can  be  approximated  by 


=  +-  +  Tui)]' 


where  I  is  an  integer  under  our  control,  and 


Ujj  =  1  -  cos  B|j(0,<l))  =  1  -  sin(\}rjj)  sin(0)  sin((|))  -  cos(\|/jj)  cos((l)) .  (9) 

Here,  we  also  used  equation  (6).  We  will  use  the  element  response  rn(0,<t>)  to  an  arrival  from 
direction  0,(t>,  as  given  by  equation  (8)  and  (9),  for  all  the  calculations  herein.  The  region  of 
applicability  is  0  <  0  <  2k,  0<^<k. 
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ELEMENT  WEIGHTS 


A  multiplicative  weight  oc  is  applied  at  general  element  location  \.,y  ,z  for 
1<^<L,  l<n<N.  The  total  effect  of  weighting  and  element  response  is  therefore  given  by 

Wn/0,<l))  =  an  P^rj^(e,(t)).  (10) 

It  is  important  to  observe  that  this  total  effective  element  response  is  factorable  in  the  variables  n 
and  I ;  it  need  not  be  factorable  in  0  and  (]). 


ARRAY  AMPLITUDE  RESPONSE 


For  notational  simplicity,  the  following  definitions  will  be  used  below: 


x^  =  27Cx^/>.,  =27tyjj /^,  Zjj  =2jcZj^ /X.,  expi(x)  =  exp(ix),  (11) 


where  X  is  the  wavelength  of  the  single-frequency  arriving  planewave.  When  the  features  above 
are  incorporated  into  the  array  amplitude  response  (2),  it  can  be  rewritten  in  the  form 


N  L 

A(0, <!))  =  ^  ^l-e  0  sin (|)  -I-  y  sin 0  sin <|)  -I-  Zj,  cos <|) 

n=U=l  *■ 

-(x^  COS0J  sin(|)^  4- y^  sin0g  sinct)^  -l-Zn  cos(t)J] , 


(12) 


where  05.<t>g  is  the  steering  direction,  and  0,<t)  is  the  planewave  arrival  direction.  Finally,  when 

the  particular  form  of  weighting  in  (10)  is  employed,  amplitude  response  (12)  further  simplifies 
to  the  desired  factored  form 


N 

A(0,  <t>)  =  ^  r^^  (0,  (}>)  exp  i[  y  sin  0  sin  <j)  -I-  z^  cos  <1)  -  a  )  x 

n=l  ' 

L 

X  y  p.  expifx.  COS0  sintj)  -  b. ) , 
e=l 


(13) 


where  parameters 
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s  sin0g  sin(l)g  +  Zji  cos^^ ,  cosS^  sincj)^  .  (14) 

The  computation  of  amplitude  response  A(0,(t))  by  means  of  (13)  requires  an  N-fold  summation 
on  n  and  a  L-fold  summation  on  i ,  which  constitutes  a  total  of  N  +  L  complex  terms.  When 
these  complex  operations  are  completed,  the  magnitude-squared  operation  yielding  power 
response  |A(0,(}))|  is  taken.  The  complex  function  expi(x)  =  exp(ix)  =  cos(x)  +  i  sin(x). 

NUMERICAL  EVALUATION  OF  DOUBLE  INTEGRAL 

The  major  computation  required  in  DI  calculation  (4)  is  the  double  integral 

n  2n 

V  H I  d(|)  sin(|)  j d0  |A(0, (^f  .  (15) 

0  0 

When  the  factorized  expression  (13)  is  substituted  in  (15),'  the  result  is  expressible  in  the  form 

n  2k  L  ^ 

V  =  J  d(t)  sin(})  I  d0  expi(x^  cos0  sin({)-b^j  x 

0  0  ^=1 

(16) 
N 

X  21  ''n  HZ  0  sin  <t)  +  z^  cos  0  -  a 

n=l 

Element  response  rj(0,(())  is  given  by  (8)  and  (9). 

The  integrand  in  (16)  is  analytic  in  0  and  (j).  Furthermore,  it  has  period  27c  in  0,  and  the  0 
integral  is  over  this  complete  interval.  This  makes  the  0  integral  in  (16)  a  good  candidate  for 
numerical  evaluation  via  the  trapezoidal  rule  [2;  pp.  53-59].  On  the  other  hand,  the  (j)  integral  is 
only  over  (0,  k),  which  is  half  the  period  in  <]).  This  suggests  the  Gauss-Legendre  rule  for 
numerical  integration  on  ([).  This  mixed  procedure  has  been  utilized  for  the  approximate 
evaluation  of  (16). 

The  specific  formula  for  M-point  Gauss-Legendre  integration  of  function  f(x)  over 
interval  (a,b)  has  the  general  form 
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(17) 


The  quantities  {tmlare  the  M  zero  locations  of  the  M-th  order  Legendre  polynomial  Pm(x),  and 
are  symmetrically  located  in  the  (-1,1)  interval.  The  M  Gauss  weights  {wm}  are  all  positive. 
These  quantities  (tm  and  Wm  for  M  up  to  768)  are  required  as  input  to  the  numerical  programs 
provided  in  the  appendices.  Gauss-Legendre  weights  and  zero-locations  are  available  in  many 
references.^  When  combined  with  the  trapezoidal  rule  for  integration  on  0,  there  follows  for 

(16),  with  sl7c4-^7rtjjj,e^  =  27tk/K, 

M  ^ 

V  =  2  IT  ^ 

m=l  k=l^=l  I 

(18) 

N  |2 

Z“n  '■n(®k''l>n,)“PVn  *“®k  +Sn  “si’m  • 

where  the  element  response  is  now  expressible  as 

=  -'l>(^u-l-|u2-l-^u^-l-'--  +  YU^j  (19) 

with 

u  s  1  -  sin\l/jj  sin0j^  ~ 

For  trapezoidal  integration,  the  interval  (0,2rc)  was  cut  into  K  panels,  each  of  width  2k/K. 
However,  the  two  usual  edge  weights  of  1/2  have  been  combined  into  a  single  unity  weight  at 
0  =  27t,  by  making  use  of  the  periodicity  in  0.  Hence,  the  summation  on  k  in  (18)  goes  from  1  to 
K,  using  unity  weights  throughout. 
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Expression  (18)  is  the  final  result  to  be  used  for  numerical  evaluation  of  the  double 
integral  for  V  in  (15),  when  the  array  is  three-dimensional.  Quantities  such  as  sinGk,  sin(l)ni,  etc., 
that  can  be  done  once  and  for  all,  have  been  precomputed  and  stored  for  access  in  (18)  as 
required,  in  order  to  minimize  the  execution  time.  A  program  for  the  evaluation  of  (18)  is 
available  in  Appendix  B  under  the  title  volumetric-grid.f. 

For  a  nonisotropic  noise  field  which  depends  only  on  polar  angle  (|),  say  with  power 
dependence  D  (([)),  the  only  changes  required  are  to  insert  D((t))  after  sin  t])  in  (15),  and  to  insert 
D((})ni)  after  sin  (l)ni  in  (18).  Of  course,  the  integral  of  product  sin((t))  D((t))  over  (0,7i)  must  also  be 
evaluated. 

SPECIALIZATION  TO  PLANAR  ARRAY,  ARBITRARY  GRID 


When  all  the  y-coordinates  {yn}  are  zero  (or  constant),  the  array  reduces  to  a  planar  array 
parallel  to  the  xz-plane;  however,  the  element  locations  can  still  be  irregular,  namely  at  each 
intersection  x.,z  forl<^<L,  l<n<N.  Every  x-coordinate  x.  is  arbitrary,  and  every 
z-coordinate  Zn  is  arbitrary  in  this  grid  structure. 

When  the  boresight  of  every  individual  element  is  perpendicular  to  the  array  plane,  angle 
Yn  in  figure  2  is  equal  to  7c/2  for  all  n.  This  makes  u  =  1  -  sin0j^  sin(|)j^  in  (20),  which  is 
independent  of  n;  therefore,  fnC^k.^n,)  (19)  is  also  independent  of  n,  allowing  it  to  be  factored 
out  of  the  n  summation  in  (18). 


Also,  when  all  the  {yn}  are  zero,  the  latter  expi  term  in  (18)  becomes  independent  of  k, 
allowing  the  entire  magnitude-squared  quantity  in  the  second  line  of  (18)  to  be  factored  out  of  the 
k  summation.  The  end  result  is  the  simplified  form 


M 


V  =  |X%sin(l)^ 


m=l 


N 

S  anexpi(znC0S(l)^  -aj 


h=l 


expi  (x^  cos 


(21) 


where  now 
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(22) 


-2|  u  +  ^u^  +  +...  +  !□* 


with 


u  =  1-sin 9j^ sin (j) .  (23) 

It  has  been  found  that  these  manipulations  and  simplifications  result  in  approximately  30 
percent  savings  in  computation  time  relative  to  the  more  general  volumetric  form  (18).  A 
program  for  simplified  form  (21)  is  listed  in  Appendix  C  under  the  title  planar-xz-grid.f,  where 
the  xz  qualifier  emphasizes  the  fact  that  the  planar  array  is  parallel  to  the  xz-plane. 


SPECIALIZATION  TO  PLANAR  ARRAY,  EQUAL  SPACINGS 


A  useful  additional  simplification  is  possible  when  the  planar  array  has  equally  spaced 
elements  in  both  the  x-  and  z-coordinates.  Specifically,  si;;:pose  that  the  element  locations  are 
given  by 


for  1  <  ^  <  L ,  Zn  =  n  d^  for  1  <  n  ^  N .  (24) 


Substitution  of  these  values  into  (21),  and  reference  to  (1 1)  and  (14),  results  in  the  modified  form 


M 


V  =  ^  V  w  sin<|> 
2  m 


m=l 


Xanexpir^d^(cos(|)n^-cos(l)3)ni 
n=l  ^  ^ 


K  ^  L  ^ 

S  K®k'1>ni|  E  Pie>!pi(x<i^(cos9^sin4i_^-cos9jSin(|>pd 

k — 1  -^—1  I 


(25) 


where  (22)  and  (23)  are  still  relevant. 

The  reduction  in  computation  time  is  now  obtained  by  taking  advantage  of  the  recursive 
capability  of  the  expi  function; 

expi[ny]  =  exp[inY]  =  exp[iY]  exp[i(n  -  1)y]  =  expi[Y]expi[(n  -  1)y]  .  (26) 
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Thus,  one  complex  multiplication  gives  the  next  term  necessary  in  the  n  or  ^  summations  in 
(25),  without  the  need  for  any  storage. 

It  has  been  found  that  these  manipulations  and  simplifications  result  in  approximately  60 
percent  additional  savings  in  computation  time  relative  to  the  more  general  planar  form  (21).  A 
program  for  simplified  planar  form  (25)  is  listed  in  Appendix  D  under  the  title  planar-xz-equal.f, 
where  the  xz-qualifier  again  emphasizes  the  fact  that  the  planar  array  is  parallel  to  the  xz-plane. 


RESULTS 


Directivity  calculations,  using  a  low-order  quadrature  formula,  are  given  in  Ref.  [4]  for  a 
9-element  line  array,  a  60-element  planar  array  and  a  21 -element  cylindrical  array.  For  arrays 
with  relatively  few  elements,  such  as  these,  there  is  little  need  to  be  concerned  with  optimal 
integration  procedures.  Computationally  efficient  integration  routines  can  provide  order-of- 
magnitude  speed  up  times,  but  this  is  of  little  importance  if  the  original  calculation  only  takes  a 
few  seconds.  However,  the  arrays  considered  here  are  large,  with  dense  element  spacings,  and 
the  frequency  and  steering  angular  dependence  of  the  directivity  must  be  examined  in  detail. 
Therefore,  it  is  essential  to  have  an  efficient  means  to  accurately  evaluate  the  two-fold  integral  of 
equation  (4). 

Table  1  below  summarizes  the  planar  and  conformal  array  lay-outs  used  for  preliminary 
numerical  directivity  calculations.  Directivity  was  calculated  at  a  number  of  frequencies  and 
assumed  a  sound  speed  of  4,900  ft/sec. 

Table  1.  Planar  and  Volumetric  Array  Input  Parameters  for  Directivity  Calculations 


Parameter 

Input  Field 

planar-xz 

equal.f 

volumetric- 

grid.f 

No.  of  Elements  along  x-axis 

L  =  200 

V 

- ^ - 

No.  of  Elements  along  z-axis 

N  =  50 

V 

X 

Element  spacing  along  x-axis 

Ax  =  3.5  in. 

V 

V 

Element  spacing  along  z-axis 

Az  =  3.6  in. 

V 

X 

Element  spacing  (arc  length) 

ds  =  3.6  in. 

X 

V 

0s-steering  angle  (azimuthal) 

0s  =  k/2  rads 

r  V 

- ^ - 

(})s-steering  angle  (polar) 

(t)s  =  7t/2  rads 

V 

— ^ — 

0-angle  increment 

A0  =  0.0001  rads 

V 

V 

(|)-angle  increment 

A(t)  =  0.0001  rads 

V 

V 

X  =  Input  not  required. 


For  both  examples,  the  total  number  of  array  elements  was  10,000.  The  angular 
increment  parameters  shown  are  necessary  in  order  to  obtain  a  reasonable  estimate  of  the 
maximum  array  response,  |A(0jjj,<t)j^)|  .  In  general,  the  number  of  angular  samples  is  related  to 
the  array  length,  and  frequency,  i.e.,  we  must  require  angular  sampling  of  the  order 
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(27) 


(A0,A4,)<^, 

J. 

where  X  is  the  acoustic  wavelength  and  .  For  the  numerical  programs  provided 

here,  this  9  and  (j)  sampling  requirement  corresponds  to  the  resection  K  >  Inlu^JX  for 
9-sampling,  where  K  is  the  program  parameter  kc,  and  M  >  -y  for  ^-sampling,  where  M 
is  the  program  parameter  me.  Additional  input  parameters,  common  to  each  program,  are  given 
in  Appendix  B. 


The  computed  directivity  is  given  in  table  2.  As  expected,  the  planar  array  geometry 
yields  slightly  greater  directivity  than  that  of  the  conformal  array;  the  differences  are  due  to  array 
curvature.  With  factorization,  it  took  154  seconds  to  compute  the  directivity  value  38.57  dB  (at 
4,000  Hz)  using  the  program  volumetric-grid.f  on  a  SUN  SPARC  Station  10.  Without 
factorization,  the  time  required  to  compute  this  value  was  on  the  order  of  1  hour,  10  minutes. 


Table  2.  Directivity  DI  (dB)  Calculated  From  planar-xz-equalf  and  volumetric-grid.f 


Frequency  (Hz) 

DI  planar-xz-equal.f  (dB) 

DI  volumetric-grid.f  (dB) 

4,000 

38.71 

38.57 

5,000 

40.63 

40.49 

6,000 

42.21 

42.06 

7,000 

43.54 

43.39 

8,000 

44.69 

44.54 

As  previously  discussed,  array  gains  in  nonisotropic  noise  fields  (having  variations  in 
polar  angle  only)  can  be  evaluated  using  the  programs  supplied  here.  In  many  ocean 
environments,  such  as  in  the  Bering  or  Norwegian  Seas,  a  distinctive  and  well-defined  notch  or 
dipole-like  directivity  pattern  is  present  in  ambient  noise.  This  variation  has  been  described  as  a 
vertical  noise  profile^  and  may  be  modeled  using  the  continuous  expression 


D(<t))  =  l-A^exp^ 


;■ 


(28) 
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Figure  3  illustrates  a  typical  nonisotropic  noise  power  dependence  for  the  parameters 
Ad  =  0.9,  which  controls  the  notch  depth,  (t>d  =  7t/2,  the  notch  location,  and  Gd  =  0.1,  the  notch 
width. 


Table  3  provides  the  array  gain  (AG)  calculated  from  planar-xz-equaLf  and  volumetric- 
grid.f  assuming  the  above  vertical  noise  power  distribution.  The  gain  over  isotropic  noise  varies 
from  6.8  dB  to  7.9  dB;  certainly,  the  vertical  noise  profile  has  a  significant  effect  on  AG 
calculations. 


Table  3.  Array  Gain  (AG)  for  Nonisotropic  Noise  Calculated  From 
the  planar-xz-equal.f  and  volumetric-grid.f  Programs 


Frequency  (Hz) 

DI  planar-xz-equal.f  (dB) 

DI  volumetric-grid.f  (dB) 

4,000 

45.54 

45.32 

5,000 

47.87 

47.65 

6,000 

49.73 

49.52 

7,000 

51.28 

51.07 

8,000 

52.60 

52.38 

APPROXIMATION  FOR  DIRECTIVITY  INDEX 


The  approximation  10  log(47tAA^)  for  the  directivity  index  DI  of  an  equi-spaced  baffled 
planar  array  is  only  applicable  for  an  unweighted  array  with  half-omni  directionality  and  an 
isotropic  noise  field.  However,  the  restriction  to  an  unweighted  array  can  be  eliminated  if  the 
area  A  is  replaced  by  the  product  of  effective  lengths  in  the  x-  and  z-directions.  That  is,  we  have, 
more  generally,  for  a  planar  array  with  spacings  dx  and  dz. 


/"47tE^  E  ^ 
DI  =  101og - ^ 

V  ^ 


where  effective  lengths  Ex  and  Ez  are  given,  respectively,  by 


(29) 
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(30) 


^  L  > 

2 

r  N 

I“n 

U=i  J 

_  p  =  H  - 

'<n=l  J 

IP? 

e=i 

’  z  z 

I“? 

n=l 

For  the  isotropic  noise  example  above,  at  4,000  Hz,  but  now  with  Hanning  weighting  in 
the  x-direction  and  flat  weighting  in  the  z-direction,  the  DI  as  found  from  the  approximation  is 
33.88  dB,  whereas  the  exact  value  turns  out  to  be  33.85  dB,  a  discrepancy  of  only  0.03  dB.  And 
for  Hanning  weighting  in  both  directions,  the  approximate  DI  is  32.12  dB,  whereas  the  exact 
value  is  32.00  dB,  a  difference  of  0.12  dB. 


POLAR  ANGLE  ^ 


Figure  3.  Typical  Variation  of  Vertical  Noise  Directionality  Assuming  Noise  Power 
Dependence  D  {^)  With  A^  =  0.9,  <1)^  =  and  =  0.1 
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APPENDIX  A  —  APPROXIMATION  TO  THE  ELEMENT  RESPONSE 


The  power  response  of  a  baffled  directive  receiving  element  can  often  be  fairly  well 
approximated  by  cos(x)  for  |x|  <  7t/2,  where  x  is  the  angle  between  the  element  boresight  angle 
and  the  planewave  arrival  angle.  For  |x|  >  7c/2,  the  response  is  ideally  zero.  The  discontinuity  in 
slope  of  the  response  at  x  =  ±7t/2  creates  a  problem  for  any  numerical  integration  technique 
which  does  not  specifically  account  for  the  locations  of  these  two  singularities. 

In  this  appendix,  we  will  derive  an  analytic  approximation  to  this  truncated  cosine,  where 
the  degree  of  fit  can  be  as  tight  as  desired.  Actually,  since  in  practice,  the  truncated  cos(x) 
response  itself  is  not  exact,  an  approximation  to  it  is  quite  satisfactory  and  acceptable.  It  will 
immediately  be  obvious  how  the  method  extends  to  a  truncated  version  of  cos'^(x),  where  power 
V  is  arbitrary. 

We  begin  by  observing  that 

cos(x)  =  exp[ln  cos(x)]  for  |x|  <  7c/2 .  (A-1) 


Alternatively, 


cos(x)  =  exp[ln(l  -  {1  -  cos(x)})]  =  exp[ln(l  -  u)] ,  (A-2) 


where 


u  s  1  -  cos(x) . 


(A-3) 


But,  from  the  expansion 


ln(l 


-u)  =  -^u  + ju^ +-|u^ +---j  for|u|<l. 


(A-4) 


that  is,  |x|  <  Jt/2,  we  immediately  obtain 
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for  all  X  in  the  extended  range  [-Tt,^],  with  u  =  1  -  cos(x).  For  large  I,  the  right-hand  side  of  (A-6) 
is  a  good  approximation  to  cos(x)  for  |x|  <  7t/2.  And,  for  |x|  >  nil,  where  u  >  1,  the  u  series 
would  tend  to  as  I  -» <»,  meaning  that  the  right-hand  side  of  (A-6)  is  approximating  zero,  as 
desired.  Furthermore,  the  right-hand  side  of  (A-6)  is  always  positive  for  all  x,  it  has  period  2ji  in 
X,  and  it  is  smooth  for  all  x.  A  sample  plot  of  the  left-hand  side  and  right-hand  side  of  (A-6)  for 
I  =  5  is  shown  in  figure  A-1.  The  transition  in  behavior  near  x  =  ±7c/2  is  now  smooth. 

The  approximation  in  (A-6)  uses  cos(x)  directly,  namely,  in  the  form  u  =  1  -  cos(x).  But, 
cos(x)  can  be  interpreted  as  the  direction  cosine  between  the  element  boresight  angle  and  the 
angle  of  arrival.  Thus,  any  other  element  response,  that  can  be  written  directly  in  terms  of  this 
quantity,  can  thereby  use  (A-6). 

For  example,  for  a  power  of  cos(x),  there  follows  from  (A-6), 

T{cos^(x) }  =  exp^-v^u  +  ^u^  +  ^u^  +  •  •  •  +  ju^  j  (A-7) 

for  |x|  e  [-7C,7c],  with  u  =  1  -  cos(x). 


A-2 


Figure  A- 1.  Sample  Plot  of  Equation  A-6 


A-3/A-4 
Reverse  Blank 


APPENDIX  B  —  Program  volumetric-grid.f 


This  program  allows  for  arbitrary  element  locations  {x^}  for  1  <  ^  <  L  and  {zn},  {yn}  for 
1  <  n  <  N .  However,  the  particular  example  programmed  utilizes  equal  spacing  dx  in  the 
x-direction  and  equal  arc  lengths  ds  in  the  polar  angle  <]),  for  simplicity.  Both  of  these 
assumptions  can  easily  be  generalized  to  fit  the  users  cases.  Then,  the  two  lines  assigning  values 
to  dx  and  ds  would  be  eliminated. 

Coding  for  Hamming  array  weighting  has  not  been  included  in  this  and  the  following 
appendix,  since  these  array  programs  can  utilize  unequal  spacings  in  the  x-  and/or  z-directions. 

It  is  certainly  possible  to  include  nonuniform  weighting,  provided  that  the  desired  shading 
function  is  properly  projected  onto  the  actual  unequal  element  locations,  accounting  for  their 
density  in  space. 


COMMON  INPUT  PARAMETERS 


K  =  400,  number  of  azimuthal  0  samples  over  (0,2k) 
M  =  192,  number  of  polar  ({)  samples  over  (0,k) 
f  =  4,000  Hz 
c  =  4,900  ft/sec 

6s  =  Jt/2,  (|)s  =  kI2\  broadside  steering 

Ax  =  3.5  inches,  Az  =  3.6  inches  for  planar  array 

Ax  =  3.5  inches.  As  =  3.6  inches  for  volumetric  array 


PROGRAM  SYMBOL 


MATH  SYMBOL  AND  EXPLANATION 


Ic 

nc 

kc 

me 

ng 

wx 

wz 

phis 

thetas 

dphi 

dtheta 

j 


L,  number  of  elements  in  x-coordinate 
N,  number  of  elements  in  z-coordinate 
K,  number  of  0  samples  over  (0,2k) 

M,  number  of  (])  samples  over  (0,7t) 

total  number  of  gauss  locations  and  weights 
storage  of  L  array  weights  for  x-coordinate 
storage  of  N  array  weights  for  z-coordinate 
(|)s,  polar  steering  angle  (radians) 

0s,  azimuthal  steering  angle  (radians) 

A(t),  phi  increment  in  search  for  maximum 
Aq,  theta  increment  in  search  for  maximum 
gauss  case  number,  1  <  j  <  12 
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implicit  reai*8 (a-h, o-z)  1  volumetric-grid. f 

real  sec, dtime, time (1 ; 2 ) 

parameter (lc=200, nc-50, kc=400, ng=2520)  1  kc  =  number  of  theta  samples 

dimension  x  (1 : Ic) , y  (1 :nc) , z (1 ;nc) ,psi (1 :nc) 

dimension  wx (1 : Ic) , wyz ( 1 : nc) , ax ( 1 : Ic) , ayz (1 :nc) 

dimension  xx (1 :ng) , ww (1 :ng) , c (1 :kc) , s (1 : kc) , cn (1 :nc) , sn (1 ; nc) 


dimension  mm(0 : 12) /O, 16, 24, 32, 48, 

1  format (d2 8 . 18)  1  case  1234 

2  format (2d25 . 16) 

3  format ("sec", 4el5 . 5) 
pi=4 . dO*atan (1 .dO) 

s  e  c=dt ime ( t ime )  1 

read  1,  xx(l:ng)  I 

read  1,  ww(l:ng)  I 

s  e  c=dt ime ( t ime ) 
print  3,  sec,  time 
print  2 

ampd= . 9d0  1 

phid=pi*.5d0  I 

sigd=.ldO  1 

freq=4000,d0  1 

speed=4900 .d0*12 .dO  1 

dx=3 . 5d0  1 

ds=3 . 6d0  1 

rc=200.d0  1 

a^O.dO  1 

b=pi  1 

phis=pi*.5d0  I 

thetas=pi* . 5d0  1 

dphi=.0001d0  I 

dtheta=.0001d0  1 

jl=8  1 

j2=10  ! 

waveiength=speed/freq 
cx=2 . dO*pi/waveiength*dx 
cyz==2  .dO*pi/wavelength*rc 

wl=0 .dO 
do  n=l,lc 

x(n)=cx*n  ! 

wx(n)=l.dO  1 

wl=wl+wx (n) 
end  do 

w2=0 .dO 
t-(nc+l) *.5d0 

da==ds/rc  1 

do  n=l,nc 
an=da* (t-n) 

y (n) =cyz*cos (an)  1 

z (n) =cyz*sin (an)  I 

psi  (n) =pi* . 5d0~an  1 

wyz(n)=l.dO  1 

w2=w2+wyz (n) 
end  do 

w3=wl*wl*w2*w2 

t=sin (phis) 
u=cos (thetas) *t 
do  n=l,lc 
ax (n)=x (n) *u 
end  do 


64,  96,  128,  192,256,  384,  512,7  68/ 

5  6  7  8  9  10  11  12  <—  j 


use  <  legendre-array-dp 
gauss-legendre  locations 
gauss-legendre  weights 


relative  power  dip,  .le.  1 
polar  angle  of  dip,  radians 
width  of  dip,  radians 

frecjuency,  hertz 
sound  speed,  inches/second 
X  spacing,  inches 
arc  length  spacing,  inches 
radius  of  array,  inches 
lower  limit  on  phi,  radians 
upper  limit  on  phi,  radians 
polar  steering  angle,  radians 
azimuth  steering  angle,  radians 
phi  increment  in  search,  radians 
theta  increment  in  search,  radians 
starting  gauss  case,  jl  .ge.  1 
ending  gauss  case,  j2  .le.  12 


X  element  locations  (arbitrary) 
x  element  weights  (arbitrary) 


element  angular  spacing,  radians 


y  element  locations  (arbitrary) 
z  element  locations  (arbitrary) 
element  polar  angles  in  y, z  plane 
yz  element  weights  (arbitrary) 
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1  theta  sampling  increment 


t=sin (thetas) *t 
u=cos (phis) 
do  n=l^nc 

ayz  (n)  =y  (n)  "^^t+z  (n)  *u 
end  do 

dt=2 .dO*pi/kc 
do  k=l^kc 
t=k*dt 
c (k) =cos  (t ) 
s (k) =sin  (t ) 
end  do 

do  n=l,nc 
t=psi (n) 
cn (n) =cos  (t) 
sn (n) =sin  (t ) 
end  do 

bal=(b+a) *.5d0 
ba2= (b-a) * . 5d0 
m22=0 
big=0 .dO 
t3-l.d0/3.d0 

do  j=l,12 

mc=mm(j)  I  me  =  number  of  phi  samples 

m22=m22+mm(  j-1) 

if  (j  .It.  jl)  go  to  9 

if  (j  .gt.  j2)  go  to  9 

v=0.d0 

two=0 .dO 

do  m=l,mc 

wm-ww  (m+m22 ) 

phim=bal+ba2*xx (m+m22) 

sm=sin  (phim) 

cm=cos  (phim) 

dm= (phim-phid) /sigd 

dm=l .dO-ampd*exp (-, 5d0*dm*dm)  1  noise  power  distribution 
vk=0.d0 

two=two+wm*sm*dm 


do  k=l,kc 
cs=c (k) *sm 
ss=s(k)*sm 
ar=0 .dO 
ai=0 .dO 
br-O.dO 
bi-O.dO 

do  n=l,lc 
wn=wx (n) 
t=x  (n)  *cs““ax  (n) 
ar=ar+wn*cos (t) 
ai=ai+wn*sin (t) 
end  do  [ n 

do  n=l,nc 

u=l  .dO-sn  (n)  *ss-cn  (n)  *cm  1  element  response: 

rn=exp (-u* (1 .dO+u* ( . 5d0+u* (t3+u* ( . 25d0+u* . 2d0) ) ) ) ) 
wn=wyz (n) *rn 
t=y (n) *ss+z (n) *cm-ayz (n) 
br=br+wn*cos (t) 
bi=bi+wn*sin (t) 
end  do  I n 
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10 


21 

4 


5 


9 


6 

7 

12 


ab= (ar*ar+ai*ai) *  (br*br+bi*bi) 
vk=vk+ab 

big)  go  to  10 


.It. 


if  (ab 
big=ab 
phib=phim 
thetab=k 
end  do  \ k 


v=v+wm*  sm*din*  vk 
end  do  Im 


v=v*ba2*dt  1  trapezoidal  and  gauss  rules 

two=two*ba2 

dfw=two*2 .d0*pi*w3/v 

diw=10.d0*logl0 (dfw) 

print  21,  v 

format ("double  integral”, d25 . 16) 


print  4,  dfw, diw 

format ("dfw,  diw", 2d25 . 16) 

dfb=two*2 .d0*pi*big/v 

dib=10.d0*logl0 (dfb) 

print  5,  dfb,dib 

format ("dfb,  dib", 2d25 . 16) 

s  e  c=dt ime ( t ime ) 
print  3,  sec,  time 
end  do  I j 

thetab=thetab*dt 
print  2 

print  6,  big, w3 

format ("big,  w3”,2d25.16) 

print  7,  phib,thetab 

format ("phib,  thetab”, 2d25 .16) 

phio=phib 
thetao=thetab 
do  i=-l,l 
phi=phio-i-dphi  *  i 
cm=cos (phi) 
sm=sin (phi) 
do  j=-l,l 

theta==thetao+dtheta*  j 

cs=cos (theta) *sm 

ss=sin (theta) *sm 

ar=0 .dO 

ai=0 .dO 

br=0.d0 

bi=0.d0 

do  n=l,lc 
wn=wx (n) 
t=x (n) *cs“ax  (n) 
ar=ar+wn*cos (t) 
ai=ai+wn*sin (t) 
end  do  1 n 

do  n=l,nc 

u=l .dO-sn (n) *ss-cn (n) *cm 

rn=exp (-u* ( 1 . dO+u* ( . 5d0+u* (t 3+u 

wn=wyz (n) *rn 

t=y (n) *ss+z (n) *cm-ayz (n) 

br=br+wn*cos (t) 

bi=bi+wn*sin (t ) 

end  do  1 n 


1  directivity  factor 
1  directivity  index 


I  big  =  maximum  power  response 
1  angles  after  coarse  search 


I  fine  search  for  maximum, 

1  starting  from  phib, thetab 


1  element  response: 
(.25d0+u*.2d0) ) ) ) ) 
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ab= (ar*ar+ai*ai) *  (br*br+bi*bi) 

if  (ab  .le.  big)  go  to  11 

big=ab 

phib=phi 

thetab=theta 

11  end  do  !  j 

end  do  1 i 

if  (abs (phib-phio) +abs  (thetab-thetao)  .gt.  O.dO)  go  to  12 
print  2 
print  8,  big 

8  format ("big”, 2d25 ,16) 

print  7,  phib,thetab  1  angles  after  fine  search 

dfb=two*2 .dO*pi*big/v  1  directivity  factor 

dib=10 .dO*loglO (dfb)  !  directivity  index 

print  5,  dfb, dib 

phib=phis 
thetab=thetas 
big=0 .dO 
k=0 

14  phio=phib 

thetao=thetab 
do  i=-l,l 
phi=phio+dphi * i 
if  (k  .eq.  0)  phi==phis 
cm=cos (phi) 
sm=sin (phi) 
do  j=-l,l 

theta=thetao+dtheta*  j 
if  (k  ,eq.  0)  theta=thetas 
cs=cos (theta) *sm 
ss=sin (theta) *sm 
ar=0 .dO 
ai=0 .dO 
br=0.d0 
bi=0.d0 

do  n=l,lc 
wn=wx (n) 
t=x (n) *cs-ax (n) 
ar=ar+wn*cos  (t) 
ai=ai+wn*sin (t) 
end  do  1 n 

do  n=l,nc 

u=l.d0-sn(n) *ss-cn(n) *cm  1  element  response: 

rn^exp (-u* (1 .dO+u* ( . 5d0+u* (t3+u* ( .25d0+u* .2d0) ) ) ) ) 
wn=wyz (n) *rn 
t=y (n) *ss+z (n) *cm-ayz (n) 
br=br+wn*cos (t) 
bi=bi+wn*sin (t) 
end  do  In 

ab= (ar*ar+ai*ai) *  (br*br+bi*bi) 
if  (ab  .le.  big)  go  to  13 
big=ab 
phib=phi 
thetab^theta 
if  (k  .eq.  0)  go  to  16 

13  end  do  I j 

end  do  I i 

16  k=k+l 

if  (k  .eq,  1)  go  to  14 


!  fine  search  for  maximum, 

1  starting  from  phis, thetas 
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if  (abs (phib-phio) +abs (thetab-thetao)  .gt.  O.dO)  go  to  14 
print  2 
print  8^  big 

print  If  phib^thetab  1  angles  near  steering  direction 

df s=two*2 . dO*pi*big/v  I  directivity  factor 

dis=10 .dO*loglO (dfs)  1  directivity  index 

print  15,  dfs,dis 
15  format ("dfs,  dis", 2d25 • 16) 

end 


B-6 


APPENDIX  C  —  Program  planar-xz-grid.f 


This  program  allows  for  arbitrary  element  locations  {  }  for  1  <  ^  <  L  and  {zn}  for 

1  <  n  <  N.  However,  the  particular  example  programmed  utilizes  two  equal  spacings  dx  and  dz, 
for  simplicity;  these  assumptions  can  easily  be  generalized  to  fit  the  users  cases.  Then,  the  two 
input  lines  assigning  values  to  dx  and  dz  would  be  eliminated. 
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implicit  real*8  (a-h, o-z)  !  pianar-xz-grid . f 

real  sec^  dt ime^  time (1:2) 
parameter (lc=200, nc=50, kc=400r ng=2520)  !  kc  =  number  of  theta  samples 
dimension  x (1 : Ic) , z  (1 :nc) , wx (1 : Ic) ^  wz (1 : nc) , ax (1 : Ic) ^ az (1 :nc) 


dimension  xx (1 :ng) , ww (1 :ng) , c 
dimension  mm(0 : 12) /O, 16, 24, 32 

1  format (d2 8 . 18)  !  case  123 

2  format (2d25. 16) 

3  format ("sec", 4el5 .5) 
pi=4 . dO  *atan ( 1 . dO ) 

sec=dtime (time) 
read  1,  xx(l:ng) 
read  1,  ww(l:ng) 
sec=dtime (time) 
print  3,  sec,  time 
print  2 

ampd= , 9d0 
phid=pi* . 5d0 
sigd=. IdO 

freq=4000 .dO 

speed=4900 .d0^12  .dO 

dx=3 . 5d0 

dz=3 . 6d0 

a=0.d0 

b=pi 

phis=pi* .5d0 

thetas=pi* . 5d0 

dphi-.OOOldO 

dtheta=, OOOldO 

jl=8 

j2-10 

wavelength=speed/ f req 
cx=2 .dO*pi/wavelength*dx 
cz=2 .dO*pi/wavelength*dz 

wl=0.d0 
do  n=l,lc 
X (n) =cx*n 
wx (n) =1 .dO 
wl=wl+wx (n) 
end  do 

w2-0.d0 
do  n=l,nc 
z (n) =cz*n 
wz  (n)  ==1  .dO 
w2=w2+wz (n) 
end  do 

w3=wl *wl *w2  *w2 

u=cos (thetas) *sin (phis) 

do  n=l,lc 

ax (n) =x (n) *u 

end  do 

u=cos (phis) 

do  n=l,nc 

az  (n)  ==z  (n)  *u 

end  do 

dt=2 .dO*pi/kc 
do  k=l,kc 
t=k*dt 


(l:kc),s(l:kc) 

.48, 64, 96, 128, 192,256,384,512,768/ 

4  5  6  7  8  9  10  11  12  <—  j 


use  <  legendre-array-dp 
gauss-legendre  locations 
gauss“legendre  weights 


relative  power  dip,  .le.  1 
polar  angle  of  dip,  radians 
width  of  dip,  radians 

frequency,  hertz 
sound  speed,  inches/second 
X  spacing,  inches 
z  spacing,  inches 
lower  limit  on  phi,  radians 
upper  limit  on  phi,  radians 
polar  steering  angle,  radians 
azimuth  steering  angle,  radians 
phi  increment  in  search,  radians 
theta  increment  in  search,  radians 
starting  gauss  case,  jl  .ge.  1 
ending  gauss  case,  j2  .le.  12 


I  X  element  locations  (arbitrary) 
1  X  element  weights  (arbitrary) 


1  z  element  locations  (arbitrary) 
I  z  element  weights  (arbitrary) 


!  theta  sampling  increment 
1  theta (k) 
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c (k) =cos  (t ) 
s (k) =sin  (t) 
end  do 

bal= (b+a) * . 5d0 
ba2=(b-a) *.5d0 
m22=0 
big=0 .dO 
t3=l.d0/3.d0 

do  j=1^12 

mc=niin(j)  1  me  =  number  of  phi  samples 

m22=m22+mm( j-1) 

if  (j  .it.  jl)  go  to  9 

if  (j  .gt.  j2)  go  to  9 

v=0.d0 

two=0 . dO 

do  m=l,mc 

wm=ww  (m+m22 ) 

phim=bal+ba2  *xx (m+m22 ) 

sm=sin  (phim) 

cm=cos  (phim) 

dm=  (phim-phid) /sigd 

dm=l .dO-ampd*exp (-. 5d0*dm*dm)  I  noise  power  distribution 
vk=0.d0 

t  wo = t  w  o +wm'*’ sm  *  dm 
br=0 .dO 
bi=0.d0 


do  n=l,nc 
wn=wz (n) 
t=z (n) *cm-az (n) 
br=br+wn*cos (t) 
bi=bi+wn*sin (t) 
end  do  I n 
bsq=br*br+bi*bi 

do  k=l,kc 
cs=c  (k) *sm 
ss=s (k) *sm 

u=l.d0-ss  1  element  response: 

rsq=exp (~2 .d0*u* (1 .dO+u* ( .5d0+u* (t3+u* ( . 25d0+u* . 2d0) ) ) ) ) 
ar=0 .dO 
ai=0 .dO 

do  n=l,lc 
wn=wx  (n) 
t=x (n) *cs-ax (n) 
ar=ar+wn*cos (t) 
ai-ai+wn*sin (t) 
end  do  In 

ab=rsq* (ar*ar+ai*ai) *bsq 
vk=vk+ab 

if  (ab  .It.  big)  go  to  10 

big=ab 

phib=phim 

thetab=k 

end  do  1 k 


v=v+wm*sm*dm*vk 
end  do  Im 

v=v*ba2*dt  I  trapezoidal  and  gauss  rules 

two=two*ba2 
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4 


5 


9 


6 

7 

12 


11 


dfw=two*2 .dO^pi^wS/v 
diw=10 . dO  ^ioglO (dfw) 
print  21/  v 

format ("double  integral"/ d25 . 16) 
print  4/  dfW/ diw 
format ("dfw,  diw", 2d25 . 16) 
dfb=two*2 .dO*pi*big/v 
dib=10.d0*logl0 (dfb) 
print  5,  dfb,dib 
format ("dfb,  dib", 2d25 .16) 

sec=dtime (time) 
print  3,  sec, time 
end  do  1 j 

thetab=thetab*dt 
print  2 

print  6,  big, w3 
format ("big,  w3",2d25.16) 
print  7,  phib,thetab 
format ( "phib,  thetab" , 2d25 . 

phio=phib  1  fine  search  for  maximum, 

thetao=thetab  I  starting  from  phib, thetab 

do  i=-l,l 

phi=phio+dphi*i 

cm^cos (phi) 

sm=sin (phi) 

do  j=-l,l 

theta=thetao+dtheta*  j 
cs=cos (theta) *sm 
ss=sin (theta) *sm 

u=l.dO-ss  1  element  response: 

rsq^exp (-2 .dO*u^ (1 .dO+u* ( . 5d0+u* (t3+u* ( . 25d0+u* . 2d0) ) ) ) ) 

ar=0 .dO 

ai=0 .dO 

br=0 .do 

bi=0 .do 

do  n=l,lc 
wn=wx (n) 
t=x (n) *cs-ax (n) 
ar=ar+wn*cos (t) 
ai=ai+wn*sin (t) 
end  do  !n 

do  n=l,nc 
wn=wz (n) 
t=z  (n) *cm-az (n) 
br=br+wn*cos (t) 
bi=bi+wn*sin (t) 
end  do  I n 

ab=rsq* (ar*ar+ai*ai) * (br*br+bi*bi) 

if  (ab  .le,  big)  go  to  11 

big=ab 

phib=phi 

thetab=theta 

end  do  1 j 

end  do  ! i 

if  (abs (phib“phio) tabs (thetab-thetao)  .gt.  O.dO)  go  to  12 

print  2 

print  8,  big 

format ("big", 2d25. 16) 

print  7,  phib, thetab  1  angles  after  fine  search 
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1  directivity  factor 
1  directivity  index 


I  big  =  maximum  power  response 
!  angles  after  coarse  search 


8 


cifb=two*2  .dO*pi*big/v 
dib=10.d0*logl0 (dfb) 
print  5,  dfb,dib 


directivity  factor 
directivity  index 


I 

I 


phib=phis  I  fine  search  for  maximum^ 

thetab=thetas  1  starting  from  phis, thetas 

big=0 .dO 
k=0 

14  phio=phib 
thetao=thetab 
do  i=-l,l 
phi=phio+dphi*i 

if  (k  .eq.  0)  phi=phis 
cm=cos (phi) 
sm=sin (phi) 
do  j=-l,l 

theta=thetao+dtheta*  j 
if  (k  .eq.  0)  theta=thetas 
cs=cos (theta) *sm 
ss=sin (theta) *sm 

u=l.d0-ss  I  element  response: 

rsq=exp (-2 .d0*u* (1 .dO+u* ( . 5d0+u* (t3+u* ( .25d0+u* .2d0) ) ) ) ) 

ar=0 .dO 

ai=0.d0 

br=0.d0 

bi=0.d0 

do  n=l, Ic 
wn=wx (n) 
t=x (n) *cs-ax (n) 
ar=ar+wn*cos (t) 
ai=ai+wn*sin  (t) 
end  do  1 n 

do  n=l,nc 
wn=wz (n) 
t=z  (n)  *cm*-az  (n) 
br=br+wn*cos (t) 
bi=bi+wn*sin (t) 
end  do  1 n 

ab=rsq* (ar*ar+ai*ai) * (br*br+bi*bi) 
if  (ab  .le.  big)  go  to  13 
big=ab 
phib=phi 
thetab=theta 
if  (k  .eq.  0)  go  to  16 
13  end  do  Ij 

end  do  1 i 
16  k-k+1 

if  (k  .eq.  1)  go  to  14 

if  (abs (phib-phio) +abs (thetab-thetao)  .gt.  O.dO)  go  to  14 
print  2 
print  8,  big 

print  1,  phib,thetab  1  angles  near  steering  direction 

dfs=two*2 .d0*pi*big/v  I  directivity  factor 

dis=10 .d0*logl0 (dfs)  I  directivity  index 

print  15,  dfs,dis 

15  format (”dfs,  dis", 2d25 . 16) 

end 


C-5/C-6 
Reverse  Blank 


APPENDIX  D  —  Program  planar«xz-equal.f 


In  this  appendix,  coding  for  Hamming  array  weighting  has  been  included  in  the  listed 
program  (as  an  option)  for  either  or  both  of  the  x-  and  z-coordinates.  However,  it  has  been 
exercised  in  only  two  of  the  numerical  examples  that  we  present  here. 


1 

2 

3 


22 

23 

c 


24 

25 
c 


implicit  real*8  (a-h, o-z)  !  planar-xz-equal . f 

real  sec, dtime, time (1 :2) 

parameter (lc=200,nc=50,kc=400,ng=2520)  1  kc  =  number  of  theta  samples 


dimension  xx (1 : ng) , ww (1 ;ng) , wx (1 : 
dimension  mm(0 : 12) /0,  16,  24,  32,  48, 
format (d2 8 . 18)  !  case  1234 
format (2d25 ,16) 
format ("sec”,  4el5 .5) 
pi=4 . dO  *atan ( 1 . dO ) 

sec=dtime (time)  1 

read  1,  xx(l:ng)  1 

read  1,  ww(l:ng)  ! 

sec=dtime (time) 
print  3,  sec, time 
print  2 

ampd= . 9d0  ! 

phid=pi*.5d0  I 

sigd=.ldO  I 

freq=4000.d0  ! 

speed=4900.d0*12.d0  I 

dx=3.5d0  I 

dz=3 . 6d0  1 

a=0,d0  I 

b=pi  1 

phis=pi*.5d0  I 

thetas^pi* . 5d0  I 

dphi=.0001d0  1 

dtheta=.0001d0  I 

jl=8  1 

j2=10  ! 

wavelength=speed/ f req 
cx=2 .dO*pi/wavelength*dx 
cz=2 .dO*pi/wavelength*dz 
bx=cx*cos (thetas) *sin (phis) 
bz=cz*cos (phis) 

wl=0.d0 

t=(lc+l) *,5d0 

if  (Ic  .gt.  1)  go  to  22 

u=0.d0 

go  to  23 

u=2 .dO*pi/ (lc-1) 

do  n=l,lc 

wx (n)=.54dO+.46dO*cos (u* (n~t) )  ! 

wx(n)=l,dO  ! 

wl=wl-l“wx  (n) 
end  do 

w2=0.d0 

t= (nc+1) * . 5d0 

if  (nc  .gt.  1)  go  to  24 

u=0 .dO 

go  to  25 

u=2 .dO*pi/ (nc“l) 

do  n=l,nc 

wz (n)=,54d0+.46d0*cos (u* (n-t) )  1 

wz(n)=l.dO  1 

w2=w2+wz (n) 
end  do 

w3=wl *wl *w2  *w2 

dt=2 .dO*pi/kc  I 


Ic) , wz (1 :nc) , c (1 :kc) , s (1 : kc) 

64,  96,  128,  192,256,  384,  512,7  68/ 

5  6  7  8  9  10  11  12  <~  j 


use  <  legendre-array-dp 
gauss-legendre  locations 
gauss-legendre  weights 


relative  power  dip,  .le.  1 
polar  angle  of  dip,  radians 
width  of  dip,  radians 

frecfuency,  hertz 
sound  speed,  inches /second 
X  spacing,  inches 
z  spacing,  inches 
lower  limit  on  phi,  radians 
upper  limit  on  phi,  radians 
polar  steering  angle,  radians 
azimuth  steering  angle,  radians 
phi  increment  in  search,  radians 
theta  increment  in  search,  radians 
starting  gauss  case,  jl  .ge.  1 
ending  gauss  case,  j2  .le.  12 


X  element  weights  (Hamming) 

X  element  weights  (arbitrary) 


z  element  weights  (Hamming) 
z  element  weights  (arbitrary) 


theta  sampling  increment 
D-2 


1  theta (k) 


do  k=l,kc 
t^k^dt 
c (k) =cos  (t) 
s (k) =sin  (t ) 
end  do 

bal=(b+a) *.5d0 
ba2=(b-a) *.5d0 
ni22=0 
big=0 .dO 
t3=l.d0/3.d0 

do  j=l, 12 

mc=niin(j)  1  inc  =  niimber  of  phi  samples 

m22=m22+mm( j“l) 

if  (j  .It.  jl)  go  to  9 

if  (j  .gt.  j2)  go  to  9 

v=0.d0 

two=0 .dO 

do  m=l^mc 
wm=ww  (m-l-m22 ) 
phim=bal+ba2*xx (m+m22) 
sm=sin  (phim) 
fx=cx*sm 

ez=cz*cos (phim) -bz 
dm=  (phim-phid) /sigd 

dm=l.d0-ampd*exp(-.5d0*dm*dm)  I  noise  power  distribution 
vk=0 . dO 

t wo=t wo +wm* sm*dm 

br=0.d0 

bi=0.d0 

r=l .dO 
q^O.dO 
ce=cos(ez) 
se=sin (ez) 
do  n=l,nc 
t=r*ce-q*se 
q=q*ce+r*se 
r=t 

wn=wz (n) 
br=br+wn*r 
bi=bi+wn*q 
end  do  I n 
bsq=br*br+bi*bi 

do  k=l,kc 

u=l  .dO-'S  (k)  *sm  1  element  response: 

rsq=exp  (-2  .dO*u*  (1  .dO+u*  ( .  5d0+u*  (t3+u*  ( .  25d0+u*  .  2d0) ) ) ) ) 
ex=fx*c (k) -bx 
ar=0 .dO 
ai=0 .dO 

r=l .dO 
q=0.d0 
ce=cos (ex) 
se=sin (ex) 
do  n=l,lc 
t=r*ce-q*se 
q=q*ce+r*se 
r=t 

wn=wx (n) 
ar=ar+wn*r 
ai=ai+wn*q 
end  do  I n 
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ab=rsq* (ar*ar+ai*ai) *bsq 
vk=vk+ab 

if  (ab  .It.  big)  go  to  10 

big=ab 

phib=phim 

thetab=k 

10  end  do  Ik 

v= V + wm  *  sm  *  dm  *  vk 
end  do  Im 

v=v*ba2*dt  !  trapezoidal  and  gauss  rules 

two=two*ba2 

dfw=two*2 .d0*pi*w3/v 

diw=10 .d0*logl0 (dfw) 

print  21,  v 

21  format ("double  integral", d2 5 .16) 

print  4,  dfw, diw 

4  format ("dfw,  diw" , 2d25 . 16) 

dfb=two*2 .d0*pi*big/v  1  directivity  factor 

dib=10 .d0*logl0 (dfb)  !  directivity  index 

print  5,  dfb,dib 

5  format ("dfb,  dib", 2d25 . 16) 

sec=dtime (time) 
print  3,  sec, time 
9  end  do  I j 

thetab=thetab*dt 
print  2 

print  6,  big, w3 

6  format ("big,  w3",2d25.16) 

print  7,  phib,thetab 

7  format ("phib,  thetab", 2d25 . 16) 

12  phio=phib 

thetao=thetab 
do  i=“l,l 
phi=phio+dphi*i 
sm=sin (phi) 
fx=cx*sm 

ez=cz*cos (phi) -bz 
br=0.d0 
bi=0.d0 

r=l.d0 
q=0.d0 
ce=cos (ez) 
seisin (ez) 
do  n=l,nc 
t=r*ce“q*se 
q=q'*'ce-Hr*se 
r=t 

wn=wz (n) 
br=br+wn*r 
bi=bi+wn*q 
end  do  1 n 
bsq=br*br+bi*bi 

do  j=-l,l 

theta=thetao+dtheta* j 

u=l . dO-sin (theta) *sm  1  element  response: 

rsq=exp (-2 .d0*u* (1 .dO+u* ( . 5d0+u* (t3+u* ( .25d0+u* .2d0) ) ) ) ) 
ex=fx*cos (theta) -bx 
ar=0 .dO 


1  fine  search  for  maximiam, 

1  starting  from  phib, thetab 


I  big  ”  maximum  power  response 
!  angles  after  coarse  search 
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ai=0 . dO 


r^l.dO 
q=0,dO 
ce=cos (ex) 
se=sin (ex) 
do  n=lrlc 
t=r *ce-q*se 
q=q*ce+r*se 
r=t 

wn=wx (n) 
ar=ar+wn*r 
ai=ai+wn*q 
end  do  In 

ab=rsq* (ar*ar+ai*ai) *bsq 
if  (ab  .le.  big)  go  to  11 
big=ab 
phib=phi 
thetab=theta 
end  do  I j 
end  do  1 i 

if  (abs  (phib-phio) tabs  (thetab-thetao)  .gt.  O.dO)  go  to  12 
print  2 
print  8,  big 
8  format ("big”, 2d25 .16) 

print  1,  phib,thetab  I  angles  after  fine  search 

dfb=two*2 .dO*pi*big/v  !  directivity  factor 

dib=10.d0*logl0 (dfb)  I  directivity  index 

print  5,  dfb,dib 

phib=phis 
thetab=thetas 
big=0 .dO 
k=0 

14  phio=phib 

thetao=thetab 
do  i=-l,l 
phi=phio+dphi*i 
if  (k  .eq.  0)  phi=phis 
sm=sin (phi) 
fx=cx*sm 

ez=cz*cos  (phi)-~bz 
br=0.d0 
bi=0.d0 

r=l.d0 
q=0.d0 

%  ce=cos(ez) 

se=sin (ez) 
do  n=l,nc 
t=r*ce-q*se 
q=q*ce+r*se 
r=t 

wn=wz (n) 
br=br+wn*r 
bi=bi+wn*q 
end  do  1 n 
bsq=br*br+bi*bi 

do  j=-l,l 

theta=thetao+dtheta*  j 
if  (k  .eq.  0)  theta=thetas 
u=l,d0-sin (theta) *sm  !  element  response: 

rsq=exp (-2 . d0*u* (1 .dO+u*  ( . SdO+u* (t3+u* ( .25d0+u* . 2d0) ) ) ) ) 
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!  fine  search  for  maximum, 

1  starting  from  phis, thetas 


1 

11 


ex=fx*cos (theta) -bx 
ar=0 .dO 
ai=0 .dO 

r=l.dO 
q=0.d0 
ce=cos (ex) 
se=sin (ex) 
do  Ic 
t=r*ce-“q*se 
q=q*ce+r*se 
r=t 

wn=wx (n) 
ar^ar+wn*r 
ai=ai'fwn*q 
end  do  ! n 

ab=rsq* (ar*ar+ai*ai) *bsq  « 

if  (ab  .le.  big)  go  to  13 
big=ab 
phib=phi 
thetab=theta 
if  (k  .eq.  0)  go  to  16 
13  end  do  Ij 

end  do  1 i 
16  k=k+l 

if  (k  .eq.  1)  go  to  14 

if  (abs (phib-phio) tabs (thetab-thetao)  .gt.  O.dO)  go  to  14 
print  2 
print  8^  big 

print  1,  phib,thetab  1  angles  near  steering  direction 

dfs=two*2 .d0*pi*big/v  1  directivity  factor 

dis=10 .d0*logl0 (df s)  I  directivity  index 

print  15,  dfs,dis 
15  format ("dfs,  dis",  2d25 . 16) 

end 
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